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Abstract 

In this short note we investigate the numerical performance of the method of artificial diffusion for 
second-order fully nonlinear Hamilton-Jacobi-Bellman equations. The method was proposed in (M. Jensen 
and I. Smears, arXiv:l 111.5423 1; where a framework of finite element methods for Hamilton-Jacobi-Bellman 
equations was studied theoretically. The numerical examples in this note study how the artificial diffusion 
is activated in regions of degeneracy, the effect of a locally selected diffusion parameter on the observed 
numerical dissipation and the solution of second-order fully nonlinear equations on irregular geometries. 



1 Introduction 



Hamilton-Jacobi-Bellman (HJB) equations, which describe the value function in the theory of optimal control, 
are fully nonlinear partial differential equations, which are of second-order if the underlying control problem 
is stochastic. One challenge arising in the numerical solution of these equations is the presence of spurious 
generalised solutions of the PDE which do not coincide with the value function. While these spurious solutions 
often possess the same regularity as the value function, they violate monotonicity properties exhibited by the 
value function. These properties lead to the concept of viscosity solutions. 

Regarding the numerical solution of HJB equations, we would like to highlight three approaches within the 
finite element methodology, which have been employed to ensure convergence to the value function. For a 
review of discrete Markov chain approximations before application of the dynamic programming principle we 
refer to 1 5 1 . For the method of vanishing moments we point to 1 3 1 . For the approach by Barles and Souganidis 
we cite the original source 1 1 1 and the recent adaption |4 | by the authors to a finite element framework. A more 
comprehensive outline of the literature is in . 
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2 Numerical Method 



Let D be a bounded Lipschitz domain in U d , d > 2. Let A be a compact metric space and 

A -> C(f2) x C(n,U d ) x C(Q) x C(Q), a — (a^kV",^") 
be continuous. Consider the bounded linear operators of non-negative characteristic form 

L a : ^(On^P-I^Q), w~-a a \w + b a -Vw + c a w 
where a belongs to A. Then 

sup||(a ,b ,c ,d ) Hc(n)xC(n,R rf )xC(n)xC(n) +SU P II Hc 2 (n)^c(n) < °° - ^ 

aeA ctE-A 

We assume that the final-time boundary data vj e C(Q) is non-negative: vj > on Q. For smooth w let i/z^ := 
sup a {L a w - d a ), where the supremum is applied pointwise. The HJB equation considered is 

-v t + Hv = in(0,T)x£3, (2a) 
v = g on (0, T) x dn, (2b) 
v = vj on {T} x Q. (2c) 

Let Vi, i = 1,2,..., be a sequence of piecewise linear shape-regular finite element spaces, whose underlying 
meshes are strictly acute. Let yf,£ = l,..., dim V,, denote the nodes of the mesh with associated hat functions 

0f. Set$ :=0f/||0f|| l 1 (tj) • Tne mesh size is denoted (Ax);. The set of time steps is Si := {s £ : k = 77 fe,}. 

Let the ^th entry of diw{s\,-) be (djH> (sf,-))^> = {w{s* +1 , y e { ) - w[s^,yp)/hi. 

For each a and i find an approximate splitting L a ~ Ef + If into linear operators 

Ef: w^> -df Aw + bf -Vw + cf w, 
If: w>-> -df Aw + bf -Vw+tf w, 

of the form a a = df + df , b a = bf +bf , c a = cf+cf &ndd a = df. To impose monotonicity, select the artificial dif- 
fusion parameters vf' e and vf' e as in [4] such that df{y^) > max{a"(yf),v"' } and df (yf) > max{a"(yf),v"' }. 

Define, for w e H l (£1), ^ e {1, . . . , N = dim V?}, 

(Efw) e :=df{/ i ){Vw,V4> £ i ) + {bf-Vw+ cf w.ft), (3a) 
(If w) e := df (yf ) (V w, V</>f > + <fcf • V w + Bf w, $\ > , (3b) 
(Cf), :=<<*?,#>. (3c) 

Obtain the numerical solution Vi(T,-) e V; by interpolation of y^. Then i>i(s!v) £ V/ at time is defined, 
inductively, by interpolating the boundary data and by 

-dt vi {sf, •) + sup(E? Vi (sf +1 , •) + \fvi (sf , •) - Cf) = 0, (4) 

a 

where the supremum is understood to be applied component-wise to the collection of vectors 

{Efvdsf +1 r) + \fvdsf,-)-Cf:aeA}. 
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Figure 1: The plot shows a locally- adapted choice of artificial diffusion over the triangular domain. The peak in 
artificial diffusion at the centre of the domain corresponds to the degeneracy of the differential operator at the 
origin. 



3 Selection of the Artificial Diffusion Parameter 

That the diffusion coefficients a" and a" in j3j are placed outside of the scalar product originates from the 
non-divergence form of the linear operators L a of HJB operators 0). This structure makes it straightforward to 

_ rv P — (X P 

implement local artificial diffusion parameters v. ' and v. ' , that is to implement a dependence on the nodal 
position yf . Indeed a change of v?' corresponds in the assembly of E" to a scalar multiplication of the ^th row 
of the stiffness matrix, see f3) . 

Changes in the optimal artificial diffusion coefficients arise from variations of the local mesh quality and the 
local mesh Peclet number. The parameters v.' and v. ' maybe chosen by studying (35) in |4| or by examining 
the assembled matrices of the unstabilised operators. This is in particular simple for the E", since only the 
signs of off-diagonal entries need be corrected to impose a local monotonicity property, thus leading to an 
algorithm which can be executed row-by- row. 

Notice that in contrast, for problems in divergence form, the stabilised diffusion coefficients a" and a" need 
to be determined in the whole domain and not just at nodal positions. 



4 An Exact Solution and Convergence Rates 

We consider a triangular spatial domain £1 with the vertices (0,-1), (v/3/2, 1/2) and (-v3/2, 1/2) and the HJB 
equation 



1 / x 2 + y 2 1 1 , , 1 Jx 2 + y 2 
-v t --\ — Av+ — Nv\ = -tj^. (5) 

2Vr-t+i 2 VT-f+i 2(r-r+i) 3/2 
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Inverse of mesh size 1/(A x) 



Figure 2: Uniform error at the initial time t = with constant ratio hi I (Ax) j 



To see that equation (5j is a HJB equation, note that the Euclidean norm of the gradient satisfies 

|Vv| = sup{/3-Vy : pen 2 with = 1}. 



(6) 



A calculation verifies that the function 



v{x,y, t) = exp 



x 2 + y 2 
T- t+l 



x 2 + y 2 
T- t+l 



is an exact solution of {5}, where boundary and final-time conditions are determined by interpolation of the 
exact solution v. The equation is to be solved on the time interval (0, 1). 

We consider a splitting which discretises the advection term explicitly with the minimum amount of diffusion 
needed for monotonicity. The remaining (linear) diffusion is incorporated in the implicit term, observing that 
this leads to a time-step restriction hi < (Ax),-. Figure [2] shows the supremum norm error with a constant ratio 
fa, /(Ax);, showing the uniform stability of the method in this setting. 

FigureJTjiHustrates a choice of artificial diffusion that is locally adapted to the Peclet number of a coarse mesh. 
It is seen that diffusion is only artificially introduced in a neighbourhood of the origin, which is where the 
operator becomes degenerate. Figure[3]illustrates the rates of convergence of the numerical solution v% (0, •) in 
the L 2 -, L°°- and H l -norms to the exact solution at the initial time, now using the constant ratio hi I (Ax)?. 



5 Eikonal Equation 



The steady-state limit of the time-dependent eikonal equation, 

-v t + \Vv\ = l, 
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Figure 3: Convergence rates of the numerical solution at the initial time t = in the L 2 -, L°°- and H 1 -norms 
with constant ratio hi I (Ax)?. 



equipped with homogeneous boundary and final-time conditions, measures the distance to the boundary. 
Due to |6j the eikonal equation belongs to the Hamilton- Jacobi-Bellman family. We consider the equation on a 
domain (Figure[4j top left) whose irregular shape leads to complicated curves on which v is not differentiable. 
The height of the domain is equal to one. 

Figure[4]compares locally-adapted choices of the artificial diffusion parameter with global choices, and illus- 
trates the effect of mesh refinement on the quality of the solution. To compare the quality of the numerical 
solutions we compare their L°°-norms — recalling that excessive numerical dissipation leads to a smearing out 
of extrema. The coarse grid solution, with 2858 internal nodes and with a global diffusion parameter of 0.05, 
only reaches a height of 9.56 x 10" 2 . In contrast, a computation on the same mesh with a local diffusion pa- 
rameter, with mean value 0.01 and standard deviation 0.005 but maximal value 0.05, leads to a height of 0.138. 
With one (two) steps of uniform refinement the number of elements increases by a factor of 4 (of 16) and the 
artificial diffusion decreases to an average value of 0.004 (of 0.002), giving the improved value of 0.147 (of 0.151) 
for the L°°-norm in the lower left (right) part of Figure [4j The L°°-norm of a reference solution on a very fine 
mesh is 0.153. 



6 A Second-Order Fully Nonlinear Equation 

The final example, on the same domain as in the previous section, is concerned with the fully nonlinear case 
of second-order. We examine the equation 

-v t + sup {-aAv} + \Vv\ = -v t + sup {-a\v + /3- Vv} = f, 

ae[a 0) ai] (a,p]E[a a ,ai\ xdB(0,l) 
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Figure 4: The top left plot shows the original mesh for the domain. The top right plot shows the numerical 
solution with globally-chosen diffusion and below is a plot with locally-chosen artificial diffusion. The bottom 
left plot is the numerical solution from one uniform refinement of the mesh and the bottom right plot shows 
the solutions after two refinements. 

where B{0, 1) is the unit ball in U 2 , T = 0.009, a = 0.045, «i = 0.09 and 

fix, y) = 529(sin(g(x, y)) + - sin(2 g{x, y)) + — sin(8g(x, y))J , 
g{x, y) = n 2 [x - 0.63) (y - 0.26) /0.07. 

The boundary and final-time conditions are homogeneous. As before, the advection term is discretised ex- 
plicitly, with the locally minimal diffusion needed for monotonicity. The possibly remaining (fully nonlinear) 
diffusion is placed in the implicit term. 

In this example the control a of the second-order term is maximised independently of the first-order control /3; 
in the sense that the optimal a may be determined without knowledge of /3. Furthermore, as Af takes locally 
either a positive or negative value only the controls ao and a\ are ever active in the HJB equation. This is an 
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Figure 5: The left plot shows the contours of the control a maximising the nonlinear second order term of the 
equation, the right plot shows the value function v. In the left plot, the black regions correspond to a = ag, 
whilst the white regions correspond to a = ct\. 

example of the Bang-bang principle. It is also reflected in the left plot in Figure[5j where the value of a is plotted. 
Black colouring signifies ao maximises the operator, whereas white colouring corresponds to a\. Observe that 
no intermediate grey values can appear. The plot of the control a mimics some of the features of the value 
function, which is plotted in the same figure on the right — in part because the Laplacian contains information 
about the curvature of the solution. 

At each time-step of the method, a semi-smooth Newton method |4| was used to solve the nonlinear discrete 
equation in l[4j, where each iteration of the algorithm involves solving a linear system. To study the perfor- 
mance of the algorithm, the HJB equation was solved on a sequence of successively refined meshes, with a 
constant set of tolerances. The linear systems were solved to a tolerance of 10" 10 by GMRES. The stopping 
criterion for Newton's method was a relative residual tolerance of 5 x 10~ 8 in the maximum norm and a con- 
vergence of the iterations requirement of 5 x 10~ 9 in the maximum norm. The sizes of the linear systems for 
the sequence of meshes were 674, 2858, 11759, 47693 and 192089. The respective average number of New- 
ton iterations for a single time step were 3, 3.67, 4.04, 4.22 and 4.86, with respective standard deviations 0, 0.48, 
0.19, 0.42 and 0.36. This demonstrates a weak dependence of the number of iterations needed for an individual 
time step on the system size, thus showing that the total number of linear systems to be solved for a complete 
computation depends principally on the number of time-steps. 
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